#Read in data
hcpd_data <- readr::read_csv("hcpd_n762_myelin_cog_data.csv")
#Column names of myelin measurements from 7 networks
seven_nets <- c('mean_wholebrain_T1wT2w', paste0('Yeo', 1:7, '_myelin')) #the first element is whole-brain
model_list <- lapply(1:length(seven_nets), function(i){
right_side <- '~ s(Age) + Sex + Scanner + numNavs_sum'
left_side <- seven_nets[[i]]
mf_form <- as.formula(paste0(left_side, gsub('s\\((.*)\\)', '\\1', right_side)))
form <- as.formula(paste0(left_side, right_side))
d <- model.frame(formula = mf_form, data = hcpd_data)
model_name <- ifelse(i == 1, 'fit_brm', paste0('fit_brm_', left_side))
if(!file.exists(paste0(model_name, '.rds'))){
fit_brm <- NULL
} else {
fit_brm <- readRDS(paste0(model_name, '.rds'))
}
return(fit_brm)
})
posterior_samples_list <- lapply(model_list, get_smooth_posterior, res, eps)
results_summaries_list <- lapply(posterior_samples_list, get_smooth_summaries)
for (i in 1:length(seven_nets)){
include_points <- TRUE #overlay points on smooth plots
cat(paste0('\n\n## ', seven_nets[[i]], '{.tabset}\n\n'))
cat(paste0('\n\n### Derivative\n\n'))
cat('Age where the median posterior derivative is at its max is marked.\n\n')
aderivative <- results_summaries_list[[i]]$derivative
anageatmax <- results_summaries_list[[i]]$max_deriv_age_posterior
asmooth <- results_summaries_list[[i]]$smooth
age_at_max_med_deriv <- aderivative[deriv1 == max(deriv1), Age]
#In order to put the smooth on the same scale as the data, we need to add back
#in the mean of the response variable that brms uses to compute the
#conditional effect of the smooth. The easiest way to do this is to use the
#output of the conditional effects output. We don't really care about the effect across
#age so we can get the output at an arbitrary age:
ce <- conditional_effects(model_list[[i]], effects = 'Age',
int_conditions = list('Age' = c(8)))
#get the model frame too, so we can plot the points to confirm this works.
mf <- brms:::model.frame.brmsfit(model_list[[i]])
#We want to reference the mean of the response variable, which I've set up to
#be in the list `seven_nets`. I'm naming the variable `response_mf_mean` to
#indicate that this is the mean of the response as determined by brms from the
#model frame (`mf`) that it uses to estimate the model.
response_mf_mean <- ce$Age[[seven_nets[[i]]]]
#We now need to add the response variable mean to the smooth estimates, but we
#only need to do that for the plot of the smooth---the derivative is already
#in the correct units (change in mylenation over age) and does not need to be
#adjusted. See below....
print(ggplot(aderivative, aes(x = Age, y = deriv1)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 1, fill = apal[2]) +
geom_segment(x = age_at_max_med_deriv,
xend = age_at_max_med_deriv,
y = aderivative[Age == age_at_max_med_deriv, lower],
yend = aderivative[Age == age_at_max_med_deriv, upper], color = apal[[1]]) +
geom_line() +
geom_line(aes(color = compared_to_max, group = 1), size = 2, alpha = .8) +
geom_point(x = age_at_max_med_deriv, y = aderivative[Age == age_at_max_med_deriv, deriv1], color = apal[[1]], size = 2) +
scale_color_manual(breaks = c('less_steep', 'as_steep'), values = apal[c(4,5)],
labels = c('Less steep', 'As steep'),
name = 'Region, as compared \nto maximum derivative, \nis...') +
jftheme)
cat(paste0('\n\n### Age at maximum derivative\n\n'))
bayesplot::color_scheme_set('red')
print(bayesplot::mcmc_areas(data.frame(max_age = anageatmax$max_age, Chain = rep(1:4, each = 2500)),
prob = 0.95,
prob_outer = 0.99) + jftheme)
cat(paste0('\n\n### Difference from derivative at steepest age\n\n'))
print(ggplot(aderivative, aes(x = Age, y = diff_from_max)) +
geom_ribbon(aes(ymin = diff_lower, ymax = diff_upper), alpha = 1, fill = apal[2]) +
geom_line(color = apal[[1]]) +
geom_line(aes(color = compared_to_max, group = 1), size = 2, alpha = .8) +
scale_color_manual(breaks = c('less_steep', 'as_steep'), values = apal[c(4,5)],
labels = c('Less steep', 'As steep'),
name = 'Region, as compared \nto maximum derivative, \nis...') +
jftheme)
cat(paste0('\n\n### Smooth estimate\n\n'))
mf_points <- NULL
if(include_points){
mf_points <- geom_point(data = mf, aes_string(y = seven_nets[[i]]), alpha = .15)
}
#This is where we finally plot the smooth! So it's the correct place to add
#in the mean of the response. We have to add it to both the line and the ribbon.
print(ggplot(asmooth[aderivative[, .(Age, compared_to_max)], on = 'Age'],
aes(x = Age, y = est + response_mf_mean)) +
mf_points +
geom_ribbon(aes(ymin = lower + response_mf_mean, ymax = upper + response_mf_mean), alpha = 1, fill = apal[2]) +
geom_line() +
geom_line(aes(color = compared_to_max, group = 1), size = 2, alpha = .8) +
scale_color_manual(breaks = c('less_steep', 'as_steep'), values = apal[c(4,5)],
labels = c('Less steep', 'As steep'),
name = 'Region, as compared \nto maximum derivative, \nis...') +
labs(y = 'est') +
jftheme)
}
Age where the median posterior derivative is at its max is marked.
Age where the median posterior derivative is at its max is marked.
Age where the median posterior derivative is at its max is marked.
Age where the median posterior derivative is at its max is marked.
Age where the median posterior derivative is at its max is marked.
Age where the median posterior derivative is at its max is marked.
Age where the median posterior derivative is at its max is marked.
Age where the median posterior derivative is at its max is marked.